Interference pattern in the collision of structures in the BEC dark matter model: 

comparison with fluids 
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In order to explore nonlinear effects on the distribution of matter during collisions within the 
Bose-Einstein condensate (BEC) dark matter model driven by the Schrodinger-Poisson system of 
equations, we study the head-on collision of structures and focus on the interference pattern for- 
mation in the density of matter during the collision process. We explore the possibility that the 
collision of two structures of fluid matter modeled with an ideal gas equation of state also forms 
interference patterns and found a negative result. Given that a fluid is the most common flavor of 
dark matter models, we conclude that one flngerprint of the BEC dark matter model is the pattern 
formation in the density during a collision of structures. 
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I. INTRODUCTION 

The discovery of the nature of dark matter is one of the 
most important problems nowadays, and a considerable 
amount of theoretical models have been proposed. One 
of the candidates that has received certain attention is 
a scalar field dark matter ultralight particle that is con- 
sistent with cosmological observations [H, Hj . At cosmo- 
logical scales, this and other related models have been 
widely studied, and one of the important issues to be 
understood about this candidate is its behavior at local 
scales related to nonlinear processes of collapse, structure 
formation and collision of structures. Some theoretical 
studies have pushed forward this model at local scales, 
especially those related to galactic rotation curves under 
various conditions • In the case of ultralight scalar 
fields, it has been shown that the density profiles of scalar 
field structures is not cuspy as opposed to Navarro Frenk 
White density profiles [1], and also predicts an adequate 
cutoff in the mass power spectrum of structures consis- 
tent with the abundance of small structures Finally, 
more recently the evolution of cosmological perturbations 
in Bose-Einstein condensate dark matter is already under 
study 0- 

In the nonlinear regime, some advances have been 
made related to the evolution of gravitating scalar field 
structures, both relativistic under general relativistic 
conditions Q and Newtonian in the low energy and weak 
gravitational field regime 

The Newtonian case has been found to be more ap- 
propriate to study the nonlinear collapse and interaction 
between structures, and is driven by the time-dependent 
Schrodinger-Poisson (SP) system of equations. The inter- 
pretation of this system of equations is that Schrodinger 
equation represents a Bose-Einstein condensate (BEC) 
of the scalar field particles at zero temperature in the 
mean fiel d ap proximation, through the Gross-Pitaevskii 
equation jl3l |. This is the reason why in the low energy 
and weak gravitational field limits, the scalar field dark 



matter model is called BEC dark matter. 

The solution of the time-dependent SP system requires 
the application of numerical methods to visualize the evo- 
lution of BEC structures and estimate any potential ob- 
servable implication. A helpful feature of the SP system 
is that it has equilibrium configurations, that is, there are 
spherically symmetric stationary solutions based on the 
assumption that the wave function depends harmonically 
in time, which in turn implies that both the gravitational 
potential and the energy density are time-independent 
• Such equilibrium configurations have been stud- 
ied dynamically and found to be stable under spherical 
perturbations and are also attractor solutions in time for 
initial perturbations with arbitrary profile ; they are 
also stable against nonspherical perturbations and 
they serve as stable structures to study how BEC dark 
matter structures behave in nonlinear situations. 

In this paper, we explore the BEC dark matter model 
du ring the head-on collision of structures as depicted first 
in |12| , where it was found that when the binary system is 
initially bounded, the resulting system is a single object, 
whereas for unbounded systems a sort of solitonic behav- 
ior was shown to happen. In both cases, it was found that 
an interference pattern was formed during the collision 
which potentially would provide predictions on the be- 
havior of baryonic matter whose stream would be driven 
in part by the dark matter gravitational potential during 
the collision of two structures. 

On the other hand, the most studied models of dark 
matter assume it can be a dust fiuid, which is an ap- 
proximate model of dark matter candidates like WIMPs. 
Thus, in this work we present a comparison of the head- 
on collision of two structures assuming on the one hand 
that the structures correspond to two equilibrium con- 
figurations of the Schrodinger-Poisson system and on the 
other hand the structures correspond to two relaxed balls 
of ideal gas. 

The goal is to determine whether or not the collision 
of two balls of gas representing two structures of WIMPs 
can also present an interference pattern during the head- 
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on collision, and thus conclude that the BEC dark matter 
has peculiar fingerprint interference patterns during the 
collision of structures that should predict observations 
that determine its viability, or determine the BEC dark 
matter should be ruled out. 

This paper is organized as follows: In Sec. II we intro- 
duce the systems of equations we solve, and the numerical 
methods we use in order to carry out the head-on colli- 
sion of two structures, both of BEC dark matter and of 
an ideal gas; in Sec. Ill we present the results of our 
simulations and compare the behavior in the two cases; 
and finally in Sec. IV we draw some conclusions. 



II. SYSTEMS OF EQUATIONS AND 
NUMERICAL METHODS 

A. The Schrodinger-Poisson system 

The system of equations. The SP system of equations 
consists of the Schrodinger equation for a wave function 
■0, with a potential that is solution of Poisson equation 
sourced by the density of probability Since we deal 
with a head-on collision, we write down the SP system 
in cylindrical coordinates as follows: 

dHJ_ ^ (2) 

where = ip{x,z.,t) and U = U{x,z,t) are the wave 
function and the gravitational potential respectively; x, z 
are the radial and axial cylindrical coordinates respec- 
tively. The third order term in Eq. ([T]) is related to a 
self-interacting term, in which A corresponds to the s- 
wave scattering length in the Gross-Pitaevskii approxi- 
mation for Bose condensates [l^, which we set to zero 
in this work, because it is irrelevant in the interference 
pattern formation studied here. This term instead was 
shown to play the role of determining the compactness 
of an equilibrium configuration [l^. Eqs. ([T]-[2]) use the 
units and scaling H = c = 1 with x — > mx, z — >■ mz, 
t — > mt and the wave function ip — > -\/47rGV', where m is 
the mass of the ultralight boson. A consequence of this 
change of units is that the mass of a system will be in 
units of [M] — M^i/m, and thus m determines the scale 
of the configurations we start with. In fact these units 
define the Kaup mass, which in the relativistic counter- 
part of the Schrodinger-Poisson system (the boson star 
case ruled by the Einstein-Klein-Gordon system of equa- 
tions) indicates the threshold between stable and unsta- 
ble bounded configurations It is worth mentioning 
that in the case of the Schrodinger-Poisson system, that 
is, the low energy and weak field limits of boson stars, the 
equilibrium configurations are all stable, unless a nega- 
tive self-interaction factor A is introduced UM. 



The evolution of the system. We consider the system 
([T][2]) to be a constrained evolution system, that is the 
Schrodinger equation is an evolution equation that satis- 
fies the Poisson equation which is a constraint. Explicitly, 
we integrate the Schrodinger equation using a finite dif- 
ferences approximation along the spatial directions and 
a method of lines for the integration in time that uses a 
third order Rungc-Kutta algorithm. We solve the con- 
straint (Poisson equation) at the intermediate steps of 
the time integrator, which provides a full coupling of the 
evolution constrained system, as shown in our conver- 
gence tests. The flow of the solution is as follows: at 
every intermediate iteration of the Runge-Kutta integra- 
tor: i) evolve the system using the Schrodinger equa- 
tion to evolve ip, ii) use such updated value of the wave 
function to source and solve the Poisson equation, iii) 
then obtain a new potential U and repeat. Usually the 
Schrodinger equation is integrated in time using unitary 
operators [l^, which is related to fully implicit meth- 
ods, however, even though we use explicit methods we 
verified the evolution is unitary as shown in [lol - [T2j . In 
order to avoid the singularity at a; = in Eqs. ([TJH) with 
our finite differences approximation, we stagger the axis 
such that we avoid the origin, and in order to achieve the 
second order convergence at the axis we use the identity 
= and code the later expression, which is a 

derivative with respect to x'^. 

Poisson equation. Eq. ([2)) is an elliptic equation for 
U which we solve using the 2D five-point stencil for the 
derivatives and a successive over-relaxation iterative al- 
gorithm with optimal acceleration parameter p^ . In 
order to impose boundary conditions we made sure the 
boundaries are far enough for the number of particles 
represented by the integral of the density of probability 
N = J \tp\'^d^x to be the same along the boundary of the 
domain and used a monopolar term of the gravitational 
field; that is, we used the value U = —N/r along the 
boundaries with r = \J x^ + z"^. At the axis we demand 
the gravitational potential to be symmetric with respect 
to the axis. 

Boundary conditions during the evolution. We use 
a sponge in the outermost region of the domain. The 
sponge is a concept used with success in the past when 
dealing with the Schrodinger equation [TtI . [18| . This 
technique consists in adding up to the potential in the 
Schrodinger equation an imaginary part. The result is 
that in the region where this is applied there is a sink 
of particles, and therefore the density of probability in 
this region will be damped out, with which we get the 
effects of an outgoing flux boundary condition [1^. The 
expression we use for the sponge profile is 

V = - '-Vo {2 + tanh[{r,k ~ r,)/6] ^ tanhirJS)} , (3) 

which is a smooth version of a step function with ampli- 
tude Vb, centered at rc and width S; rjk = x^j. + zj^, is 
the radius of a given point of the discretized domain; the 
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minus sign guarantees the decay of the number of parti- 
cles at the outer parts of our integration domain, that is, 
the imaginary potential behaves as a sink of particles. 

Initial data. In order to reduce our parameter space 
to specific but illustrative cases, we only consider the 
equal mass head-on collision case. Thus we construct the 
solution of a spherically symmetric equilibrium ground 
state configuration in spherical coordinates as done in 

[13, [13 ■ III the units used the mass of each config- 
uration is iV = 2.06 with a radius containing the 95% 
of the total mass is xg^ =3.93 [l^. We obtain a wave 
function tpsph and gravitational potential Usph- We then 
interpolate such wave function at two different places of 
the two-dimensional grid (0, ±zo) centered along the z- 
axis and define a global wave function ip = "01 + "V^Zj 
where ipi and V'2 arc the wave functions obtained from 
the spherical solution located at two different places of 
the 2D grid. We choose zq such that the two configu- 
rations are far enough one from the other so that the 
interference term < ipi,ip2 > lies near to round off er- 
ror values, in order to consider the two blobs have the 
adequate phase; the magnitude of this term decreases 
exponentially with the initial separation of the blobs. In 
this way, we solve the Poisson equation sourced by the 
energy density p = I'i/'p. Then we have initial data for 
two superposed ground state equilibrium configurations 
in our axially symmetric domain. 

Analysis. In order to analyze the kinetic development 
of the collision it is important to track global quantities 
during the evolution, such as kinetic, gravitational and 
total energy. We do this by calculating the expectation 
values 



K = --J V*VVrf^a;, 



(4) 
(5) 



which are respectively the expectation values of the ki- 
netic and gravitational energies. These quantities are 
important at determining the state of the system at any 
time during the evolution of the system. For instance, the 
value of the total energy E = K + W indicates whether 
we account with a bounded system or not, and the very 
important virial theorem relation 2K + W = 0, which is 
nearly satisfied when the system gets virialized and re- 
laxed through whatever channels available, for instance, 
the emission of scalar field bursts, the so-called gravita- 
tional cooling process [l9j . 



in (20[ . This scheme allows us to evolve the position, ve- 
locity, density, pressure and internal energy of the fluid 
elements. We explore two different scenarios of ideal. In 
the first one we set the internal energy and pressure equal 
to zero, and in the second one we use the more general 
ideal gas equation of state p = (T — l)pu, where p is the 
pressure of the gas, p is its density, F is a constant and 
u is the internal energy of the gas. 

Initial data. The initial values required for our simu- 
lations are the positions, velocities, density, pressure and 
internal energy of the particles in the fluid. 

We start choosing the density proflle given by the 
Plummer model: 



for the fluid without pressure and the density profile: 



(6) 
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(7) 



for the fluid with pressure. 

We construct one of the structures using an 
acceptance-rejection method to generate the position of 
the particles in such a way that they satisfy the cor- 
responding density profile. We set the initial velocities 
of the particles equal to zero and set the parameters 
M = 2.0622 and R = 3.93. Finally, we set the value 
of the internal energy equal to a constant (u = for the 
pressureless fluid and u = 0.05M/i? in the other case) 
and recover the pressure from the equation of state. Such 
configurations are not in equilibrium, therefore we evolve 
them numerically and let them settle down into an equi- 
librium state. 

In order to generate the initial data corresponding to 
the head-on collision of two bumps of fluid, we take the 
resulting equilibrium configuration described in the pre- 
vious paragraph, make two copies of it and place them 
at a given separation along the z-axis. Then a boost on 
each of the configurations is applied along the head-on 
axis. The initial positions are z = ±10.0 and the initial 
velocity is = ±2.6875/2.0622. This value was chosen 
such that the total energy of the system is equal to zero, 
although other values of the boost were explored and the 
results are qualitatively the same. 



III. RESULTS 



B. The SPH algorithm for the fluid 

We want to study the head-on collision of two struc- 
tures of ideal gas. In order to do that, we need to solve 
Euler's equations coupled to Newtonian gravity. With 
this in mind, we implemented a smoothed particle hydro- 
dynamic (SPH) code based on the formulation described 



We perform a series of head-on collisions of two initial 
balls in both cases, the one driven by the SP system of 
equations and the fiuid driven by hydrodynamics. We set 
up configuration for the SP system and for the fluid that 
are dynamically comparable in time, size and separation. 
The reason is that there is no way to identify exactly 
an initial configuration corresponding to the quantum 
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mechanical system driven by the Schrodinger equation 
(which is dispersive) with the configurations constructed 
for an ideal gas. 

We define the total energy of the system a.s E = K+W, 
where K and W are the expectation values of the kinetic 
gravitational energies defined in (|4][5]), calculated with 
the time-dependent wave function that is being calcu- 
lated on the fly of the Schrodingcr-Poisson system, and 
E = K + U + W is the total energy, that is the sum of the 
kinetic, internal and gravitational energies in the case of 
the fluid. The total energy depends strongly on the value 
of the initial head-on momentum the initial balls have, 
so that we choose this parameter to be the one we tune 
in order to obtain our initial binary conflgurations. 

The parameter space we explore corresponds to initial 
configurations with different values of the total energy E, 
which allows to explore situations in which the system is 
both bounded {E < 0) and unbounded {E > 0). 

Wc consider there is a pattern formation if there is a 
sequence of strips of the density of matter of the size of 
the order of the initial radii of the initial configurations 
as shown to happen during the collision of two blobs of 
BEC [13. 

A. Bose condensate initial configurations 

We choose two ground state spherically symmetric 
equilibrium configurations that are originally virialized, 
that is, they satisfy separately the condition 2K + W = 
and superpose them into the 2d grid (for details in the 
construction of equilibrium configurations see [l8|). We 
apply momentum along the head-on direction in the fol- 
lowing way: assume ijj = + "02 is the wave function 
resulting from the superposition of the two equilibrium 
configurations, where ipi and '02 represent the wave func- 
tion of each of the two equilibrium configurations cen- 
tered along the head-on axis z at ±zo, being ipi and ■02 
centered at — zq and zq respectively. Thus we add up 
momentum to each of the configurations defining a new 
wave function ■0 = e^^'i/ji + e~'^"''02- After this process 
we solve Poisson equation and start the evolution. 

According to the Gross-Pitaevskii theory, BECs be- 
have like classical waves, and as such interference phe- 
nomena are expected. Experiments showing interfer- 
ence patterns are well known [2l|. In these experiments 
two separate Bose condensate configurations are released 
from their respective traps and allow them to evolve 
freely; eventually the two densities of probability inter- 
fere and produce a total density of the form 

IV'l+V'21' = |^l|' + |V2|' + 2 < V^1,V^2 > COS , 

.(8) 

where m is the mass parameter in Schrodinger's equation, 
and d is the initial distance between the two condensates, 
and $ is the initial phase difference between the two con- 
densated given by ^0 = ipi+e'^^ip2, which is nonzero when 



there is a nonzero head-on initial momentum. The den- 
sity of probability ([5]) predicts a spacing between two 
consecutive fringes of constructive interference given by 



which in our units corresponds to A = 2TTt/d = 27r</(2zo). 
This type of pattern formation has also been verified 
with simulations of laboratory systems (see e.g. [2^). 
We have to remind the reader that an important differ- 
ence between the conditions in these experiments and our 
system is that we actually never release the condensates 
from their respective traps which are the gravitational 
potentials the density of probability generates, instead, 
we allow the traps to interact with each other. Neverthe- 
less we use the formula to measure the distance between 
fringes and find that for big values of the initial head-on 
momentum p^, that is, a big value of $ in the experi- 
ments, the law is nearly valid, whereas for small phases 
it does not hold. We interpret this result as a conse- 
quence of the fact that for the cases where the head-on 
momentum is small, a merger is actually expected to oc- 
cur and the interaction between the traps becomes impor- 
tant, whereas in the high momentum a nearly solitonic 
trespassing effect happens. 

The parameter space we consider is shown in Table 
U and snapshots of pattern formation are shown in Fig. 
[T] Notice that the width of the energy density in the 
patterns is even thicker than the initial size of the equi- 
librium configurations at initial time. 
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3.99 


2.3 


1.5 
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~ 6.5 


2.01 


1.8 


2.5 


4- 


~ 4.3 


1.27 


1.2 



TABLE I: Parameters of the initial configurations for the 
BEC case. We also show the value of the distance between 
fringes of interference both, as predicted by formula (j9]) for 
laboratory experiments and the one we measure. 



We perform our calculations using the domain x € 
[0,20], z e [—20,20] and a uniformly discretized grid, 
with constant time resolution At and the same resolu- 
tion along both spatial directions Ax = Az = Axz 
and a sponge radius with Tc = 17 and 6 = 1. In or- 
der to validate our numerical results we show in Fig. [2] 
a criterion that shows our results converge. As men- 
tioned above, we use a method of lines with second 
order stencils along the spatial directions and a third 
order accurate integrator in time, so that we expect 
the code to converge to second order at least. Since 
we have no exact solution to our problem to compare 
with, or a constraint of the PDF system of equations 
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FIG. 1: We show the best snapshot of the density \ip'^ along 
the z axis showing the interference pattern. We also show the 
initial value of the density, in order to compare the size of the 
interference pattern. 



wc should expect to be satisfied, we ean only practice a 
self-convergence test, that is, we use three different res- 
olutions Aszi = 0.2, IS.XZ2 = 0.13333 = Axzi/1.5 and 
Axz^ = 0.08888 = Aa;z2/1.5 and track the maximum of 
the energy density pbec = IV'Pj which corresponds to the 
infinity norm of pf,ec because it is always non-negative. 
Given the ratio between our resolutions is 1.5, we expect 
the maximum of the density to satisfy 

max{phf,c[using Axi]) — max{pbec[using Aa;2]) ^ 
max{phec[using Aa:2]) — max{phec[using Asa]) 

(10) 

where Q is the order of convergence expected to hold, 
which is at least two in our case because such is the or- 
der of accuracy of our stencils along the spatial directions 
and at most three, which is the accuracy of our time inte- 
grator. In Fig. [5] we also show the value of Q. In all our 
calculations we kept the factor At/ Axz^ = 0.1 constant 
for all the resolutions, which maintains the evolution sta- 
ble, accurate and convergent. 

B. Fluid configurations 

We evolve the initial configurations of fluid integrating 
Euler's equations using a predictor-corrector scheme. We 
compute the total, kinetic, potential and internal energy 
of the system in order to monitor the behavior of our 
numerical evolution. 

The first case we analyzed is the pressureless fiuid. In 
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FIG. 2: (Top) In order to show the convergence of our algo- 
rithms we track the maximum of the energy density during 
the evolution for the particular case = 1.5, which would 
be equivalent to estimate the infinity norm of the energy den- 
sity. In our calculations we used three different resolutions 
Axi = 0.2, Ax2 = Aa;i/1.5 and Axs = Aa;2/1.5. (Bottom) 
We show the convergence factor, that is 1.5*^, which in the- 
ory should be 1.5^ — 2.25 for a second order convergent im- 
plementation. We show how much our calculations approach 
second order convergence; the convergence factor oscillates 
usually near or around the theoretical value due to phases in 
the scalar (the maximum of the density in our case) for the 
various resolutions. A big peak also appears approximately 
by the time the density approaches its maximum, and we see 
how convergence is lost after around t ~ 11 and decreases 
clearly afterwards by the time the blobs reach the sponge re- 
gion, where the calculations are not expected to converge since 
there Schrodinger's equation has been modified with an arti- 
ficial sink of particles, and also the unitarity of the evolution 
in the numerical domain is lost. 



Fig. [3]we show snapshots at different times of the density 
as function of the z-position of the particles. The time 
of collision is around t ~ 7.2 and it is clear that no inter- 
ference pattern is generated during such a collision. In 
Fig. |3]we show the kinetic, potential, internal and total 
energy of the system, in which it can be noticed that no 
internal energy is being involved in the model. 

The second case is the ideal gas case with pressure. As 
before, in Fig. [5] we present the density as function of 
the position for different times and in Fig. [Slthe energies 
of the system. In this scenario, again, we were unable to 
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track any interference patterns. These figures illustrate 
the role played by the pressure, which produces the den- 
sity profiles to be less steep compared to the pressureless 



case. 
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FIG. 5: In this plot we show snapshots at different times of 
the density as function of the z coordinate of the structures 
of fluid. As in the pressureless case, the evolution shows the 



FIG. 3: In this plot we show snapshots at different times of 
the density as function of the z coordinate of the structures 
of an ideal gas with zero pressure. The collision takes place 
after t ~ 7.2. The evolution shows that the fluids cross each 
other without showinp' a.nv nattern of interference. 




b3 



-10 - 





1 1 1 










Kinetic 








■ ■ ■ Potential 








Internal 




5 




— - Total 










Energy 
o 








-5 


I.I 







I , I I I , I , I 

5 10 15 20 



t 

FIG. 6: In this plot we show the total, kinetic, potential and 
thermal energy for the head-on collision of fluid with pressure. 
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FIG. 4: In this plot we show the total, kinetic, potential and 
thermal energy for the head-on collision of the pressureless 
fluid conflgurations. The peak appears approximately by the 
time the density of the fluid reaches its maximum. 

We have also performed simulations with various val- 
ues of the initial linear momentum in the head-on direc- 
tion and thus various values of the total energy of the 
system. Wc have also carried out our simulations with 
various numbers of particles and in all the cases the re- 



sults are qualitatively the same. 

IV. CONCLUSIONS 

Wc have solved the SP system as the equations describ- 
ing the evolution of gravitating Bose condensates, which 
represents the BEC dark matter model and its properties 
at local scales. We have focused on the head-on collision 
of two equilibrium ground state configurations, and stud- 
ied the interference pattern formation of the density of 
the configuration during the collision. 
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In order to investigate whether or not other type of flu- 
ids may show a similar interference pattern, we also stud- 
ied the collision of two spherical configurations made of 
an ideal gas fluid, both with and without pressure using 
SPH techniques for different values of the total energy. 

We found that the pattern formation during the col- 
lision of structures does not happen for the fluid. This 
leads to the conclusion that a fingerprint of the BEC dark 
matter model is the presence of interference patterns dur- 
ing the collision of structures. That is, if evidence of such 
interference patterns is not found, the BEC dark matter 



model should be ruled out. 
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